In [47]:
import sys
sys.setrecursionlimit(10000)
import warnings
warnings.filterwarnings('ignore', category=DeprecationWarning)

import os
os.environ['GNUMPY_IMPLICIT_CONVERSION'] = 'allow'
print os.environ.get('GNUMPY_IMPLICIT_CONVERSION')
allow

In [29]:
import cPickle
import gzip

from breze.learn.data import one_hot
from breze.learn.base import cast_array_to_local_type
from breze.learn.utils import tile_raster_images

import climin.stops
import climin.initialize
from climin import mathadapt as ma

from breze.learn import hvi
from breze.learn.hvi import HmcViModel
from breze.learn.hvi.energies import (NormalGaussKinEnergyMixin, DiagGaussKinEnergyMixin, MlpDiagGaussKinEnergyMixin)
from breze.learn.hvi.inversemodels import MlpGaussInvModelMixin, MlpInvAcceptProbModelMixin

from matplotlib import pyplot as plt
from matplotlib import cm

import numpy as np
from scipy.stats import norm as normal_distribution

#import fasttsne

from IPython.html import widgets
%matplotlib inline

import theano
theano.config.compute_test_value = 'ignore'#'raise'
from theano import (tensor as T, clone)
In [4]:
datafile = '../mnist.pkl.gz'
# Load data.                                                                                                   

with gzip.open(datafile,'rb') as f:                                                                        
    train_set, val_set, test_set = cPickle.load(f)                                                       

X, Z = train_set                                                                                               
VX, VZ = val_set
TX, TZ = test_set

Z = one_hot(Z, 10)
VZ = one_hot(VZ, 10)
TZ = one_hot(TZ, 10)

X_no_bin = X
VX_no_bin = VX
TX_no_bin = TX

# binarize the MNIST data
np.random.seed(0)
VX = np.random.binomial(1, np.tile(VX, (5, 1))) * 1.0
TX = np.random.binomial(1, np.tile(TX, (5, 1))) * 1.0
X  = np.random.binomial(1, X) * 1.0

print VX.shape

image_dims = 28, 28

X_np, Z_np, VX_np, VZ_np, TX_np, TZ_np, X_no_bin_np, VX_no_bin_np, TX_no_bin_np = X, Z, VX, VZ, TX, TZ, X_no_bin, VX_no_bin, TX_no_bin
X, Z, VX, VZ, TX, TZ, X_no_bin, VX_no_bin, TX_no_bin = [cast_array_to_local_type(i) 
                                                        for i in (X, Z, VX,VZ, TX, TZ, X_no_bin, VX_no_bin, TX_no_bin)]
print X.shape
(50000L, 784L)
(50000L, 784L)

\\srv-file.brml.tum.de\nthome\cwolf\code\2015-christopherwolf-msc\breze\learn\base.py:39: UserWarning: Implicilty converting numpy.ndarray to gnumpy.garray
  warnings.warn('Implicilty converting numpy.ndarray to gnumpy.garray')

In [5]:
fig, ax = plt.subplots(figsize=(9, 9))

img = tile_raster_images(X_np[:64], image_dims, (8, 8), (1, 1))
ax.imshow(img, cmap=cm.binary)
Out[5]:
<matplotlib.image.AxesImage at 0x2476fcf8>
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [6]:
fast_dropout = False

if fast_dropout:
    class MyHmcViModel(HmcViModel, 
                   hvi.FastDropoutMlpBernoulliVisibleVAEMixin, 
                   hvi.FastDropoutMlpGaussLatentVAEMixin, 
                   MlpDiagGaussKinEnergyMixin,
                   MlpGaussInvModelMixin,
                   MlpInvAcceptProbModelMixin):
        pass

    kwargs = {
        'p_dropout_inpt': .1,
        'p_dropout_hiddens': [.2, .2],
    }

    print 'yeah'

else:
    class MyHmcViModel(HmcViModel, 
                   hvi.MlpBernoulliVisibleVAEMixin, 
                   hvi.MlpGaussLatentVAEMixin, 
                   DiagGaussKinEnergyMixin,  # MlpDiagGaussKinEnergyMixin, DiagGaussKinEnergyMixin
                   MlpGaussInvModelMixin):  # MlpInvAcceptProbModelMixin
        pass
    kwargs = {}


batch_size = 500
#optimizer = 'rmsprop', {'step_rate': 1e-4, 'momentum': 0.95, 'decay': .95, 'offset': 1e-6}
#optimizer = 'adam', {'step_rate': .5, 'momentum': 0.9, 'decay': .95, 'offset': 1e-6}
optimizer = 'adam', {'step_rate': 0.001}

# This is the number of random variables NOT the size of 
# the sufficient statistics for the random variables.
n_latents = 20
n_hidden = 200

m = MyHmcViModel(X.shape[1], n_latents, 
                 [n_hidden, n_hidden], ['softplus'] * 2, 
                 [n_hidden, n_hidden], ['rectifier'] * 2, 
                 [n_hidden, n_hidden], ['rectifier'] * 2,
                 n_hmc_steps=3, n_lf_steps=4,
                 n_z_samples=1,
          optimizer=optimizer, batch_size=batch_size, allow_partial_velocity_update=False, perform_acceptance_step=False,
          compute_transition_densities=True, consider_aux_inv_model_inpt_constant=False,
          **kwargs)

climin.initialize.randomize_normal(m.parameters.data, 0.1, 1e-1)
#m.parameters.__setitem__(m.hmc_sampler.step_size_param, 0.2)
#m.parameters.__setitem__(m.kin_energy.mlp.layers[-1].bias, 1)
\\srv-file.brml.tum.de\nthome\cwolf\code\ml_support\theano\theano\scan_module\scan_perform_ext.py:135: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility
  from scan_perform.scan_perform import *

Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [7]:
old_best_params = None
#print m.score(TX)
print m.parameters.data.shape
(616685,)

In [8]:
FILENAME = 'hvi_gen2_recog2_aux2_late20_hid200_pretr_3hmc_04lf.pkl'

# In[5]:
#old_best_params = None
f = open(FILENAME, 'rb')
np_array = cPickle.load(f)
old_best_params = cast_array_to_local_type(np_array)
f.close()
print old_best_params.shape
(616685L,)

In [9]:
m.parameters.data = old_best_params.copy()
old_best_loss = m.score(VX)
print old_best_loss
compiled score function
garray(92.49337768554688)

Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [10]:
#print m.score(VX)
print m.score(TX)
print m.score(X[:50000])
garray(91.78533935546875)
garray(91.50416564941406)

In [13]:
print m.parameters.view(m.init_recog.mlp.layers[2].bias)
garray([ -2.69359316e-05,   2.17193499e-01,   2.57873267e-01,
          2.54865848e-02,   3.21485698e-02,   1.42345145e-01,
          3.25903833e-01,   3.86436284e-03,  -3.07249665e-01,
          1.31776467e-01,   4.22098301e-02,   4.19317596e-02,
         -2.16690749e-01,   1.17295846e-01,   1.00071228e-03,
          1.31734535e-01,   4.02299732e-01,  -7.67451897e-02,
         -3.49123150e-01,  -3.65785271e-01,  -5.12402654e-01,
         -2.67633110e-01,  -8.97625089e-02,  -5.54986954e-01,
         -7.60390639e-01,  -2.47353390e-01,  -2.16502279e-01,
         -8.21632445e-01,  -2.28132725e-01,  -1.61989808e-01,
         -7.29750693e-02,  -1.58693790e-01,  -8.00007582e-02,
         -3.07665110e-01,  -3.12527031e-01,  -1.13484986e-01,
         -4.54903662e-01,  -3.02022338e-01,  -4.16113734e-01,
         -3.55050474e-01])

In [None]:
#m.parameters.__setitem__(m.hmc_sampler.step_size_param, m.parameters.view(m.hmc_sampler.step_size_param) * np.sqrt(4.0/2.0))
#m.parameters.__setitem__(m.init_recog.mlp.layers[2].bias, cast_array_to_local_type(np.array([ 0.14096579,  0.53084856, -1.75648117, -2.5653615 ])))
#m.parameters.__setitem__(m.kin_energy.variance_parameter, cast_array_to_local_type(-0.7*np.ones_like(m.kin_energy.variance_parameter)))
In [11]:
print 0.1 * m.parameters.view(m.hmc_sampler.step_size_param) ** 2 + 1e-8
garray([ 0.097611])

In [47]:
print m.parameters.view(m.hmc_sampler.partial_velocity_parameter)
print np.tanh(m.parameters.view(m.hmc_sampler.partial_velocity_parameter))
[ 0.48037896]
[ 0.44654706]

In [16]:
#print m.estimate_nll(TX, 100)   # 93.5946457865
#print m.estimate_nll(TX, 1000)  # 92.4632834801
print m.estimate_nll(TX, 10000)  # 91.6807054881
#print m.estimate_nll(TX, 100)
#print m.estimate_nll(TX, 100)
91.6807054881

In [None]:
m._f_loss, m._f_dloss = m._make_loss_functions()
In [None]:
if False:
    full_f_dloss = m._f_dloss
    only_aux_inv_mask = np.zeros_like(m.parameters.data)
    for i in range(3):
        auxiliary_inv_weight_param_range = m.parameters._var_to_slice[m.auxiliary_inv_model.mlp.layers[i].weights]
        auxiliary_inv_bias_param_range = m.parameters._var_to_slice[m.auxiliary_inv_model.mlp.layers[i].bias]
        print auxiliary_inv_weight_param_range, auxiliary_inv_bias_param_range
        only_aux_inv_mask[slice(*auxiliary_inv_weight_param_range)] = 1.0
        only_aux_inv_mask[slice(*auxiliary_inv_bias_param_range)] = 1.0

    def only_aux_inv_loss(params, data, *args):
        return cast_array_to_local_type(only_aux_inv_mask) * full_f_dloss(params, data, *args)

    m._f_dloss = only_aux_inv_loss
In [None]:
from theano.printing import debugprint
#debugprint(m._f_dloss.theano_func)
In [None]:
print m.parameters._var_to_shape
print
print m.parameters._var_to_slice
In [None]:
grad_array = np.zeros((100, len(m.parameters.data)))
for i in range(100):
    grad_array[i, :] = m._f_dloss(m.parameters.data, X[:625])
In [None]:
mean_grad = grad_array.mean(axis=0)
print abs(mean_grad).mean()
print mean_grad.min(), mean_grad.argmin()
print mean_grad.max(), mean_grad.argmax()

# 0.0068310414746
# -8.18554486275 196327
# 14.5694446664 197230
In [None]:
x = (abs(mean_grad) > 5).astype('int32')
matchings_indices = [ i for i, value in enumerate(x) if value == 1 ]
print matchings_indices
# [196327, 197226, 197230, 197270, 197450, 197566, 197750, 197862, 197866, 197938]
In [None]:
std_grad = grad_array.std(axis=0)
print std_grad.mean()
print std_grad.min(), std_grad.argmin()
print std_grad.max(), std_grad.argmax()
print mean_grad[std_grad.argmax()]

# 0.0173799498656
# 0.0 2
# 179.490609337 197230
# 14.5694446664
In [None]:
print (abs(mean_grad) > 1).sum()
print (abs(mean_grad) > 2).sum()
print (abs(mean_grad) > 10).sum()
# 310
# 55
# 2
In [None]:
consconst_mean_grad = mean_grad
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [12]:
TARGET_FILENAME = 'hvi_gen2_recog2_aux2_late2_hid200_pretr_3hmc_04lf_accept_with_prob'
FILETYPE_EXTENSION = '.pkl'
old_best_params = None

m.optimizer = 'adam', {'step_rate': 0.0001}

max_passes = 200
max_iter = max_passes * X.shape[0] / batch_size
n_report = X.shape[0] / batch_size

stop = climin.stops.AfterNIterations(max_iter)
pause = climin.stops.ModuloNIterations(n_report)

# theano.config.optimizer = 'fast_compile'

for i, info in enumerate(m.powerfit((X_no_bin,), (VX,), stop, pause, eval_train_loss=False)):
    print i, m.score(X[:10000]).astype('float32'), info['val_loss'], np.exp(m.parameters.view(m.kin_energy.variance_parameter).as_numpy_array()), 0.1*m.parameters.view(m.hmc_sampler.step_size_param).as_numpy_array()**2 + 1e-8
    if i == 0 and old_best_params is not None:
        if info['best_loss'] > old_best_loss:
            info['best_loss'] = old_best_loss
            info['best_pars'] = old_best_params
    
    if info['best_loss'] == info['val_loss']:
        f = open(TARGET_FILENAME + FILETYPE_EXTENSION, 'wb')
        cPickle.dump(m.parameters.data, f, protocol=cPickle.HIGHEST_PROTOCOL)
        f.close()
Compiled loss functions
0 125.787994385 129.204925537 [[ 0.41872048  0.54157557]] [ 0.00429188]
1 125.700050354 129.02116394 [[ 0.41884893  0.54148064]] [ 0.00427786]
2 125.635520935 128.99005127 [[ 0.41893382  0.5408776 ]] [ 0.00430395]
3 125.702125549 129.097640991 [[ 0.41876827  0.5407894 ]] [ 0.00432126]
4 125.669219971 129.05368042 [[ 0.41890845  0.54116857]] [ 0.00428372]
5 125.711868286 129.033889771 [[ 0.41859346  0.54130981]] [ 0.00429574]
6 125.680221558 129.070037842 [[ 0.41845459  0.54124944]] [ 0.00431583]
7 125.769996643 129.068161011 [[ 0.41851547  0.54113257]] [ 0.00431906]
8 125.715049744 129.073669434 [[ 0.41923245  0.54066699]] [ 0.00428504]
9 125.710144043 129.078521729 [[ 0.4194272   0.54099841]] [ 0.00423742]
10 125.701545715 128.973892212 [[ 0.41881293  0.54096107]] [ 0.00429782]
11 125.769424438 129.05821228 [[ 0.41935051  0.54116289]] [ 0.00422406]
12 125.654418945 129.02796936 [[ 0.41894578  0.54098577]] [ 0.00427729]
13 125.651443481 129.122283936 [[ 0.41886223  0.54103566]] [ 0.00428311]
14 125.674850464 129.037414551 [[ 0.41880886  0.54082077]] [ 0.00429767]
15 125.661735535 128.925231934 [[ 0.41909034  0.54053826]] [ 0.00427935]
16 125.731269836 129.134017944 [[ 0.41916903  0.54052183]] [ 0.00427687]
17 125.624221802 129.034469604 [[ 0.41878183  0.54078389]] [ 0.00431016]
18 125.598175049 128.967086792 [[ 0.4187665   0.54139709]] [ 0.00427775]
19 125.699707031 129.098083496 [[ 0.41885474  0.54080097]] [ 0.00430295]
20 125.650047302 129.023132324 [[ 0.41875759  0.54071846]] [ 0.00431853]
21 125.6222229 128.988052368 [[ 0.41939495  0.54020999]] [ 0.00427771]
22 125.678321838 129.003540039 [[ 0.41945092  0.53978796]] [ 0.00429626]
23 125.659507751 129.040786743 [[ 0.41945785  0.53987478]] [ 0.00428932]
24 125.718833923 128.980758667 [[ 0.41967164  0.53990255]] [ 0.0042682]
25 125.730895996 129.108688354 [[ 0.4200643  0.5394477]] [ 0.00426023]
26 125.668006897 129.078292847 [[ 0.41984022  0.53917055]] [ 0.00429275]
27 125.58367157 128.961380005 [[ 0.41997398  0.53929712]] [ 0.00427142]
28 125.704437256 129.048812866 [[ 0.41995493  0.53909362]] [ 0.00428491]
29 125.63609314 128.985580444 [[ 0.41990689  0.53904487]] [ 0.00429375]
30 125.670021057 129.114883423 [[ 0.41983354  0.53966974]] [ 0.00426374]
31 125.81552124 129.033538818 [[ 0.41999365  0.53939909]] [ 0.00426541]
32 125.680168152 129.078872681 [[ 0.41976063  0.53958991]] [ 0.00427545]
33 125.704055786 129.079330444 [[ 0.41953951  0.54032653]] [ 0.0042538]
34 125.682762146 129.067489624 [[ 0.4190491   0.54005073]] [ 0.00432092]
35 125.702285767 128.989349365 [[ 0.41919986  0.54030109]] [ 0.00429396]
36 125.763771057 129.071182251 [[ 0.4190236   0.54041791]] [ 0.00429952]
37 125.662994385 129.002960205 [[ 0.41959741  0.54011279]] [ 0.00425982]
38 125.652137756 129.002304077 [[ 0.41951721  0.54022806]] [ 0.00425848]
39 125.688873291 129.091583252 [[ 0.41945367  0.54032859]] [ 0.00426086]
40 125.666923523 129.076965332 [[ 0.41930529  0.54044838]] [ 0.0042663]
41 125.688346863 129.079833984 [[ 0.41950401  0.54047766]] [ 0.00423768]
42 125.663658142 129.106384277 [[ 0.41935813  0.54061208]] [ 0.00424848]
43 125.641082764 129.078079224 [[ 0.4193962   0.54013024]] [ 0.0042766]
44 125.674758911 129.017669678 [[ 0.41956712  0.53975669]] [ 0.00428664]
45 125.674285889 128.95262146 [[ 0.41971745  0.53949857]] [ 0.00428807]
46 125.659873962 129.033035278 [[ 0.41944917  0.53955948]] [ 0.00431449]
47 125.741508484 129.051071167 [[ 0.41949778  0.53968733]] [ 0.00429515]
48 125.722259521 129.137161255 [[ 0.41943885  0.53988436]] [ 0.00429543]
49 125.725418091 129.000579834 [[ 0.41944985  0.54034689]] [ 0.00425277]
50 125.676856995 129.100143433 [[ 0.41920718  0.54018256]] [ 0.00429376]
51 125.676948547 129.068130493 [[ 0.41915961  0.54004696]] [ 0.00430387]
52 125.674598694 129.051498413 [[ 0.41932441  0.5393544 ]] [ 0.00433162]
53 125.712562561 129.030029297 [[ 0.41950053  0.54007448]] [ 0.00427063]
54 125.669456482 129.025680542 [[ 0.4194845   0.54040406]] [ 0.00425432]
55 125.703155518 129.045883179 [[ 0.41929055  0.54033452]] [ 0.00427577]
56 125.683525085 129.076187134 [[ 0.41919654  0.54047998]] [ 0.00426623]
57 125.688446045 128.93788147 [[ 0.41894301  0.54072526]] [ 0.00427641]
58 125.737571716 128.954666138 [[ 0.41865397  0.54108426]] [ 0.00427745]
59 125.782493591 129.065338135 [[ 0.4188323   0.54109687]] [ 0.00426395]
60 125.760025024 129.021652222 [[ 0.41897492  0.54045895]] [ 0.00429033]
61 125.773582458 128.984222412 [[ 0.41892341  0.54088453]] [ 0.00427257]
62 125.683143616 129.055038452 [[ 0.41865075  0.54064327]] [ 0.00431483]
63 125.683486938 129.019790649 [[ 0.4186525   0.54098116]] [ 0.00430161]
64 125.774398804 129.073318481 [[ 0.41824378  0.541481  ]] [ 0.00431203]
65 125.729545593 129.01411438 [[ 0.41847055  0.54157692]] [ 0.00427268]
66 125.704147339 129.070510864 [[ 0.4182093  0.5423955]] [ 0.00425364]
67 125.684135437 128.953460693 [[ 0.41798467  0.54206988]] [ 0.00428905]
68 125.762756348 129.032791138 [[ 0.41750667  0.54234141]] [ 0.00430811]
69 125.681144714 129.085678101 [[ 0.41757449  0.54254142]] [ 0.00429175]
70 125.707168579 129.073471069 [[ 0.41762061  0.54246048]] [ 0.00429197]
71 125.738883972 129.05241394 [[ 0.41784025  0.54211906]] [ 0.00429344]
72 125.717811584 128.988540649 [[ 0.41783519  0.54249473]] [ 0.00428278]
73 125.667671204 128.929412842 [[ 0.41776222  0.54242628]] [ 0.00428478]
74 125.684181213 129.018341064 [[ 0.41799563  0.54227201]] [ 0.00427596]
75 125.788955688 128.981170654 [[ 0.41798061  0.54208225]] [ 0.0042871]
76 125.727409363 128.928573608 [[ 0.41806946  0.54222405]] [ 0.0042658]
77 125.663619995 129.047073364 [[ 0.4183854   0.54180927]] [ 0.00426202]
78 125.701950073 128.996246338 [[ 0.4185626   0.54042367]] [ 0.00433033]
79 125.75151825 128.981643677 [[ 0.41867735  0.54049664]] [ 0.00431821]
80 125.64617157 128.961654663 [[ 0.41897759  0.54078279]] [ 0.00427085]
81 125.760231018 128.913619995 [[ 0.41866218  0.5408698 ]] [ 0.00429217]
82 125.626182556 128.910842896 [[ 0.41841191  0.54173603]] [ 0.00426133]
83 125.718170166 128.995697021 [[ 0.41779275  0.54203886]] [ 0.00430011]
84 125.715934753 128.944061279 [[ 0.41783427  0.542065  ]] [ 0.00428937]
85 125.674209595 128.980377197 [[ 0.41818877  0.54198959]] [ 0.0042632]
86 125.632659912 128.982467651 [[ 0.41823366  0.54197696]] [ 0.00426552]
87 125.661010742 128.984222412 [[ 0.41801395  0.54212846]] [ 0.0042733]
88 125.679420471 128.903244019 [[ 0.41791671  0.54227722]] [ 0.00425925]
89 125.754859924 129.060760498 [[ 0.41852321  0.54199838]] [ 0.00423311]
90 125.708572388 128.957092285 [[ 0.41859446  0.54183513]] [ 0.00423726]
91 125.745223999 128.944213867 [[ 0.41868129  0.54163771]] [ 0.0042375]
92 125.760124207 129.04524231 [[ 0.41898379  0.54133762]] [ 0.0042344]
93 125.680374146 128.971496582 [[ 0.4188706   0.54065571]] [ 0.00429083]
94 125.669418335 129.041397095 [[ 0.41865731  0.54080333]] [ 0.00429881]
95 125.760055542 129.026763916 [[ 0.41871693  0.54003566]] [ 0.00434398]
96 125.757072449 128.960830688 [[ 0.41891267  0.54030528]] [ 0.00430517]
97 125.70892334 129.015869141 [[ 0.41906199  0.54034219]] [ 0.004292]
98 125.71395874 129.067138672 [[ 0.41909044  0.54085468]] [ 0.00424892]
99 125.75037384 128.95880127 [[ 0.41919476  0.54036425]] [ 0.00426886]
100 125.698623657 129.024551392 [[ 0.41940077  0.53987706]] [ 0.00427967]
101 125.986923218 129.191223145 [[ 0.41940417  0.53999054]] [ 0.00425678]
102 125.724647522 129.001190186 [[ 0.42024201  0.53999314]] [ 0.00418223]
103 125.734672546 128.959472656 [[ 0.42005251  0.53947015]] [ 0.00423913]
104 125.731872559 128.958816528 [[ 0.42010622  0.53959319]] [ 0.00423674]
105 125.709487915 128.965637207 [[ 0.41992454  0.53917942]] [ 0.00427689]
106 125.862045288 129.063537598 [[ 0.41998456  0.53975402]] [ 0.00422668]
107 125.721336365 129.035202026 [[ 0.4193815   0.53979237]] [ 0.00427419]
108 125.724693298 128.994934082 [[ 0.41919671  0.54008852]] [ 0.00426923]
109 125.729644775 129.013397217 [[ 0.41917308  0.54004458]] [ 0.00427537]
110 125.748718262 128.938934326 [[ 0.41924846  0.540521  ]] [ 0.00424027]
111 125.724525452 129.043533325 [[ 0.41913353  0.54017396]] [ 0.00427089]
112 125.691818237 128.99937439 [[ 0.41930914  0.53991085]] [ 0.00428001]
113 125.784111023 128.963760376 [[ 0.4189914  0.5401692]] [ 0.00428203]
114 125.717262268 128.968276978 [[ 0.41931384  0.53987671]] [ 0.0042789]
115 125.702522278 128.931152344 [[ 0.41924317  0.53938867]] [ 0.00430871]
116 125.62789917 128.937133789 [[ 0.41931839  0.53891704]] [ 0.00432795]
117 125.63747406 128.973403931 [[ 0.41983882  0.53912858]] [ 0.00427694]
118 125.645187378 128.995285034 [[ 0.41973669  0.53915708]] [ 0.00427607]
119 125.753570557 129.006530762 [[ 0.41980231  0.53918369]] [ 0.00428059]
120 125.702819824 128.913116455 [[ 0.41968995  0.53898412]] [ 0.00429093]
121 125.725135803 128.907592773 [[ 0.41989746  0.53812922]] [ 0.00432323]
122 125.75869751 128.962539673 [[ 0.4199735   0.53836538]] [ 0.00429974]
123 125.765769958 129.015213013 [[ 0.42045139  0.53815058]] [ 0.00426201]
124 125.748062134 129.046691895 [[ 0.42052826  0.53756176]] [ 0.00430648]
125 125.648071289 128.989624023 [[ 0.42079689  0.537751  ]] [ 0.0042769]
126 125.754112244 129.079498291 [[ 0.42064493  0.53665861]] [ 0.00433896]
127 125.659385681 129.00541687 [[ 0.42081187  0.53746507]] [ 0.00428889]
128 125.730873108 129.006072998 [[ 0.42084949  0.53720817]] [ 0.00429186]
129 125.793235779 128.988555908 [[ 0.42094287  0.53759332]] [ 0.00426626]
130 125.738731384 128.957901001 [[ 0.42104254  0.53717446]] [ 0.00428043]
131 125.763648987 128.950653076 [[ 0.421456    0.53703731]] [ 0.0042623]
132 125.630012512 128.92388916 [[ 0.42125732  0.53753036]] [ 0.00424066]
133 125.678222656 128.964004517 [[ 0.42119447  0.53713123]] [ 0.00427011]
134 125.686073303 129.011749268 [[ 0.42085278  0.53737739]] [ 0.00428361]
135 125.717071533 129.074493408 [[ 0.42084016  0.53726159]] [ 0.00429295]
136 125.64276886 129.010818481 [[ 0.42105772  0.53760351]] [ 0.00424348]
137 125.723457336 129.089950562 [[ 0.42086269  0.53730338]] [ 0.00428478]
138 125.646095276 129.03125 [[ 0.42056942  0.53727587]] [ 0.00430859]
139 125.626350403 129.056762695 [[ 0.42035504  0.53817724]] [ 0.00426133]

\\srv-file.brml.tum.de\nthome\cwolf\code\climin\climin\util.py:150: UserWarning: Argument named f is not expected by <class 'climin.adam.Adam'>
  % (i, klass))

In [None]:
train_result_params = m.parameters.data.copy()
#m.parameters.data = info['best_pars']
#m.score(VX)
In [None]:
m.parameters.data = train_result_params
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [17]:
f_z_init_sample = m.function(['inpt'], m.init_recog.sample(), numpy_result=True)
f_z_sample = m.function(['inpt'], m.hmc_sampler.output, numpy_result=True)
f_gen = m.function([m.gen.inpt], m.gen.sample(), numpy_result=True)
f_gen_rate = m.function([m.gen.inpt], m.gen.rate, numpy_result=True)
f_joint_nll = m.function(['inpt'], m.joint_nll, numpy_result=True)
In [18]:
curr_pos = T.matrix('current_position')
curr_vel = T.matrix('current_velocity')
norm_noise = T.matrix('normal_noise')
unif_noise = T.vector('uniform_noise')

new_sampled_vel = m.hmc_sampler.kin_energy.sample(norm_noise)
updated_vel = m.hmc_sampler.partial_vel_constant * curr_vel + m.hmc_sampler.partial_vel_complement * new_sampled_vel
performed_hmc_steps = m.hmc_sampler.perform_hmc_steps(curr_pos, curr_vel)
hmc_step = m.hmc_sampler.hmc_step(curr_pos, curr_vel, np.float32(0), norm_noise, unif_noise)
lf_step_results = m.hmc_sampler.simulate_dynamics(curr_pos, curr_vel, return_full_list=True)

f_pot_en = m.function(['inpt', curr_pos], m.hmc_sampler.eval_pot_energy(curr_pos), numpy_result=True)
f_kin_en = m.function(['inpt', curr_vel], m.kin_energy.nll(curr_vel).sum(-1), numpy_result=True)
f_perform_hmc_steps = m.function(['inpt', curr_pos, curr_vel], 
                                T.concatenate([performed_hmc_steps[0], performed_hmc_steps[1]], axis=1))
f_hmc_step = m.function(['inpt', curr_pos, curr_vel, norm_noise, unif_noise], 
                        T.concatenate([hmc_step[0], hmc_step[1]],axis=1), on_unused_input='warn')
f_kin_energy_sample_from_noise = m.function(['inpt', norm_noise], new_sampled_vel)
f_updated_vel_from_noise = m.function(['inpt', curr_vel, norm_noise], updated_vel)
f_perform_lf_steps = m.function(['inpt', curr_pos, curr_vel],
                               T.concatenate([lf_step_results[0], lf_step_results[1]], axis=0))
\\srv-file.brml.tum.de\nthome\cwolf\code\2015-christopherwolf-msc\breze\arch\util.py:630: UserWarning: theano.function was asked to create a function computing outputs given certain inputs, but the provided input variable at index 5 is not part of the computational graph needed to compute the outputs: uniform_noise-for-gpu.
To make this warning into an error, you can pass the parameter on_unused_input='raise' to theano.function. To disable it completely, use on_unused_input='ignore'.
  on_unused_input=on_unused_input, updates=updates)

In [19]:
z1_by_z0_0 = T.grad(lf_step_results[0][-1, 0, 0], curr_pos)
z1_by_z0_1 = T.grad(lf_step_results[0][-1, 0, 1], curr_pos)
z1_by_z0 = T.concatenate([z1_by_z0_0, z1_by_z0_1], axis=0)
f_z1_by_z0 = m.function(['inpt', curr_pos, curr_vel], z1_by_z0, numpy_result=True, on_unused_input='warn')

zT_by_z0_0 = T.grad(performed_hmc_steps[0][-1, 0, 0], curr_pos)
zT_by_z0_1 = T.grad(performed_hmc_steps[0][-1, 0, 1], curr_pos)
zT_by_z0 = T.concatenate([zT_by_z0_0, zT_by_z0_1], axis=0)
f_zT_by_z0 = m.function(['inpt', curr_pos, curr_vel], zT_by_z0, numpy_result=True, on_unused_input='warn')
In [20]:
f_z_init_mean = m.function(['inpt'], m.init_recog.mean, numpy_result=True)
f_z_init_var = m.function(['inpt'], m.init_recog.var, numpy_result=True)

f_v_init_var = m.function(['inpt'], T.extra_ops.cpu_contiguous(m.kin_energy.var), numpy_result=True)

full_sample = m.hmc_sampler.sample_with_path()
f_full_sample = m.function(['inpt'], T.concatenate([full_sample[0], full_sample[1]], axis=1))
In [21]:
pos = T.matrix('pos')
vel = T.matrix('vel')
updated_vel = T.matrix('updated_vel')
time_step = T.vector('time_step')

aux_inpt_replacements = {m.auxiliary_inv_model_inpt['position']: pos, 
                         m.auxiliary_inv_model_inpt['time']: T.cast(time_step[0], dtype='float32')}
if 'velocity' in m.auxiliary_inv_model_inpt:  # only use the updated velocity if it is part of the model
    aux_inpt_replacements.update({m.auxiliary_inv_model_inpt['velocity']: updated_vel})

aux_inv_model_var = clone(m.auxiliary_inv_model.var, replace=aux_inpt_replacements)
aux_inv_model_mean = clone(m.auxiliary_inv_model.mean, replace=aux_inpt_replacements)
aux_inv_model_nll = clone(m.auxiliary_inv_model.nll(vel).sum(-1), replace=aux_inpt_replacements)

if 'velocity' not in m.auxiliary_inv_model_inpt:
    f_aux_inv_var = m.function(['inpt', pos, time_step], aux_inv_model_var, numpy_result=True)
    f_aux_inv_mean = m.function(['inpt', pos, time_step], aux_inv_model_mean, numpy_result=True)
    f_aux_inv_nll = m.function(['inpt', pos, time_step, vel], aux_inv_model_nll, numpy_result=True)
else:
    f_aux_inv_var = m.function(['inpt', pos, updated_vel, time_step], aux_inv_model_var, numpy_result=True)
    f_aux_inv_mean = m.function(['inpt', pos, updated_vel, time_step], aux_inv_model_mean, numpy_result=True)
    f_aux_inv_nll = m.function(['inpt', pos, updated_vel, time_step, vel], aux_inv_model_nll, numpy_result=True)
        
final_vel_inpt_replacements = {m.final_vel_model_inpt['position']: pos, 
                               m.final_vel_model_inpt['time']: T.cast(m.n_hmc_steps, dtype='float32')}

final_vel_model_var = clone(m.final_vel_model.var, replace=final_vel_inpt_replacements)
final_vel_model_mean = clone(m.final_vel_model.mean, replace=final_vel_inpt_replacements)
final_vel_model_nll = clone(m.final_vel_model.nll(vel).sum(-1), replace=final_vel_inpt_replacements)

f_v_final_var = m.function(['inpt', pos], final_vel_model_var, numpy_result=True)
f_v_final_mean = m.function(['inpt', pos], final_vel_model_mean, numpy_result=True)
f_v_final_model_nll = m.function(['inpt', pos, vel], final_vel_model_nll, numpy_result=True)

f_kin_energy_nll = m.function(['inpt'], m.kin_energy.expected_nll, numpy_result=True)
In [22]:
f_init_recog_nll = m.function(['inpt'], m.init_recog.expected_nll.sum(-1), numpy_result=True)
In [None]:
pos = T.matrix()
f_pot_en_deriv = m.function(['inpt', pos], m.hmc_sampler.eval_pot_energy_deriv(pos))
f_pot_en_2nd_deriv0 = m.function(['inpt', m.hmc_sampler.pot_energy_inpt], T.grad(m.hmc_sampler.pot_energy_deriv[0, 0], m.hmc_sampler.pot_energy_inpt))
f_pot_en_2nd_deriv1 = m.function(['inpt', m.hmc_sampler.pot_energy_inpt], T.grad(m.hmc_sampler.pot_energy_deriv[0, 1], m.hmc_sampler.pot_energy_inpt))
In [23]:
print f_init_recog_nll(VX).mean()
init_var = f_z_init_var(VX)
print init_var.mean()
print init_var.max()
print init_var.min()
print
print f_joint_nll(TX).mean()
0.843194
0.105364
0.655432
0.00406902

92.8733

Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [24]:
fig, axs = plt.subplots(2, 3, figsize=(27, 18))

### Original data

O = (X_no_bin_np[:64])[:, :784].astype('float32')
img = tile_raster_images(O, image_dims, (8, 8), (1, 1))
axs[0, 0].imshow(img, cmap=cm.binary)

O2 = (X_np[:64])[:, :784].astype('float32')
img = tile_raster_images(O2, image_dims, (8, 8), (1, 1))
axs[1, 0].imshow(img, cmap=cm.binary)

### Reconstruction

#z_sample = f_z_sample((X[:64]))
z_init_sample = cast_array_to_local_type(f_z_init_sample((X[:64])))
z_sample = f_perform_hmc_steps((X[:64]), 
                               z_init_sample, 
                               f_kin_energy_sample_from_noise((X[:64]), 
                                                              cast_array_to_local_type(np.random.normal(size=(64, m.n_latent)).astype('float32')))
                               )[-1, :64, :]

R = f_gen_rate(z_sample)[:, :784].astype('float32')
img = tile_raster_images(R, image_dims, (8, 8), (1, 1))
axs[0, 1].imshow(img, cmap=cm.binary)

Rinit = f_gen_rate(z_init_sample)[:, :784].astype('float32')
img = tile_raster_images(Rinit, image_dims, (8, 8), (1, 1))
axs[0, 2].imshow(img, cmap=cm.binary)

R2 = f_gen(z_sample)[:, :784].astype('float32')
img = tile_raster_images(R2, image_dims, (8, 8), (1, 1))
axs[1, 1].imshow(img, cmap=cm.binary)

Rinit2 = f_gen(z_init_sample)[:, :784].astype('float32')
img = tile_raster_images(Rinit2, image_dims, (8, 8), (1, 1))
axs[1, 2].imshow(img, cmap=cm.binary)
Out[24]:
<matplotlib.image.AxesImage at 0xf0cca080>
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [25]:
fig, axs = plt.subplots(1, 2, figsize=(18, 9))

prior_sample = cast_array_to_local_type(np.random.randn(64, m.n_latent))

S = f_gen_rate(prior_sample)[:, :784].astype('float32')
img = tile_raster_images(S, image_dims, (8, 8), (1, 1))
axs[0].imshow(img, cmap=cm.binary)

S2 = f_gen(prior_sample)[:, :784].astype('float32')
img = tile_raster_images(S2, image_dims, (8, 8), (1, 1))
axs[1].imshow(img, cmap=cm.binary)

#S3 = f_gen_rate(prior_sample)[:, :784].astype('float32')
#img = tile_raster_images(S, image_dims, (8, 8), (1, 1))
#axs[2, 2].imshow(img, cmap=cm.nipy_spectral)
Out[25]:
<matplotlib.image.AxesImage at 0xee7ea128>
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [23]:
dim1 = 0
dim2 = 1
In [24]:
unit_interval_positions = np.linspace(0.025, 0.975, 20)
positions = normal_distribution.ppf(unit_interval_positions)
print unit_interval_positions
print positions

latent_array = np.zeros((400, m.n_latent))
latent_array[:, dim2] = -np.repeat(positions, 20)  # because images are filled top -> bottom, left -> right (row by row)
latent_array[:, dim1] = np.tile(positions, 20)
        
fig, axs = plt.subplots(1, 1, figsize=(24, 24))

F = f_gen_rate(cast_array_to_local_type(latent_array)).astype('float32')

img = tile_raster_images(F, image_dims, (20, 20), (1, 1))
#axs.imshow(img, cmap=cm.nipy_spectral)
axs.imshow(img, cmap=cm.binary)
[ 0.025  0.075  0.125  0.175  0.225  0.275  0.325  0.375  0.425  0.475
  0.525  0.575  0.625  0.675  0.725  0.775  0.825  0.875  0.925  0.975]
[-1.95996398 -1.43953147 -1.15034938 -0.93458929 -0.75541503 -0.59776013
 -0.45376219 -0.31863936 -0.18911843 -0.06270678  0.06270678  0.18911843
  0.31863936  0.45376219  0.59776013  0.75541503  0.93458929  1.15034938
  1.43953147  1.95996398]

Out[24]:
<matplotlib.image.AxesImage at 0x119fa5710>
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [26]:
data_set = 'train'  # 'train', 'validation', 'test'

if data_set == 'train':
    L = f_z_sample(X)
    L_init = f_z_init_sample(X)
    class_vec = Z[:].argmax(1)
elif data_set == 'validation':
    L = f_z_sample(VX)
    L_init = f_z_init_sample(VX)
    class_vec = VZ[:].argmax(1)
elif data_est == 'test':
    L = f_z_sample(TX)
    L_init = f_z_init_sample(TX)
    class_vec = TZ[:].argmax(1)
else:
    print 'unknown data set, must be one of: train, validation, test'
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [31]:
dim1 = 0
dim2 = 1
In [32]:
fig, axs = plt.subplots(1, 2, figsize=(18, 9))

normal_cdf_rescale = True

if normal_cdf_rescale:
    end_points = normal_distribution.cdf(L)
    init_points = normal_distribution.cdf(L_init)
else:
    end_points = L
    init_points = L_init

axs[0].scatter(end_points[:, dim1], end_points[:, dim2], c=class_vec, lw=0, s=10, alpha=.4)
axs[1].scatter(init_points[:, dim1], init_points[:, dim2], c=class_vec, lw=0, s=10, alpha=.4)

cax = fig.add_axes([0.95, 0.2, 0.02, 0.6])
cax.scatter(np.repeat(0, 10), np.arange(10), c=np.arange(10), lw=0, s=300)
cax.set_xlim(-0.1, 0.1)
cax.set_ylim(-0.5, 9.5)
plt.yticks(np.arange(10))
plt.tick_params(axis='x', which='both', bottom='off', top='off', labelbottom='off')
cax.tick_params(axis='y', colors='white')
for tick in cax.yaxis.get_major_ticks():
    tick.label.set_fontsize(14)
    tick.label.set_color('black')
    
cax.spines['bottom'].set_color('white')
cax.spines['top'].set_color('white') 
cax.spines['right'].set_color('white')
cax.spines['left'].set_color('white')

axs[0].set_title('After HMC steps')
axs[1].set_title('Initial recognition model')

if normal_cdf_rescale:
    axs[0].set_xlim(0, 1)
    axs[0].set_ylim(0, 1)
    axs[1].set_xlim(0, 1)
    axs[1].set_ylim(0, 1)
else:
    axs[0].set_xlim(-3, 3)
    axs[0].set_ylim(-3, 3)
    axs[1].set_xlim(-3, 3)
    axs[1].set_ylim(-3, 3)
In []:
if n_latents == 2:
    if normal_cdf_rescale:
        fig, axs = plt.subplots(1, 1, figsize=(24, 24))

        F = f_gen_rate(cast_array_to_local_type(latent_array)).astype('float32')
        img = tile_raster_images(F, image_dims, (20, 20), (1, 1))
        #axs.imshow(img, cmap=cm.nipy_spectral)
        axs.imshow(img, cmap=cm.binary)
        axs.scatter(578.0*end_points[:, dim1], 578.0*(1 - end_points[:, dim2]), c=class_vec, lw=0, s=10, alpha=0.7)

        axs.set_xlim([0, 578])
        axs.set_ylim([578, 0])
    else:
        print 'Currently only available, if the scatter plot used rescaling by the normal CDF.'
        print 'Set normal_cdf_rescale=True in the previous cell'
In [33]:
fig, axs = plt.subplots(4, 5, figsize=(20, 16))
colors = cm.jet(np.linspace(0, 1, 10))
for i in range(5):
    axs[0, i].scatter(init_points[Z[:].argmax(1) == i, dim1], init_points[class_vec == i, dim2], c=colors[i], lw=0, s=5, alpha=.2)
    axs[1, i].scatter(end_points[Z[:].argmax(1) == i, dim1], end_points[class_vec == i, dim2], c=colors[i], lw=0, s=5, alpha=.2)
    axs[0, i].set_title(str(i) + ' before HMC')
    axs[1, i].set_title(str(i) + ' after HMC')
    axs[2, i].scatter(init_points[Z[:].argmax(1) == (5+i), dim1], init_points[class_vec == (5+i), dim2], c=colors[5+i], lw=0, s=5, alpha=.2)
    axs[3, i].scatter(end_points[Z[:].argmax(1) == (5+i), dim1], end_points[class_vec == (5+i), dim2], c=colors[5+i], lw=0, s=5, alpha=.2)
    axs[2, i].set_title(str(5+i) + ' before HMC')
    axs[3, i].set_title(str(5+i) + ' after HMC')
    for j in range(4):
        if normal_cdf_rescale:
            axs[j, i].set_xlim(0, 1)
            axs[j, i].set_ylim(0, 1)
        else:
            axs[j, i].set_xlim(-3, 3)
            axs[j, i].set_ylim(-3, 3)
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [37]:
X_index = 2  # index=0 -> 5, index=1 -> 0, index=2 -> 4, index=3 -> 1, index=24 -> underlined 1, index=39 -> ugly 6
from_data_set = X
from_data_set_no_bin = X_no_bin
num_repeats = 1000
np.set_printoptions(precision=4, suppress=True)

X_array = np.array([from_data_set[X_index, :]]).astype('float32')
X_array_no_bin = np.array([from_data_set_no_bin[X_index % (from_data_set_no_bin.shape[0]), :]]).astype('float32')
fig, axs = plt.subplots(1, 2, figsize=(6, 3))

img = tile_raster_images(X_array, image_dims, (1, 1), (1, 1))
axs[0].imshow(img, cmap=cm.binary)
img = tile_raster_images(X_array_no_bin, image_dims, (1, 1), (1, 1))
axs[1].imshow(img, cmap=cm.binary)
Out[37]:
<matplotlib.image.AxesImage at 0xe42beba8>
In [38]:
repeated_X = cast_array_to_local_type(np.tile(X_array, (num_repeats, 1)))

full_sample = f_full_sample(repeated_X).astype('float32')
z_samples = full_sample[:, :num_repeats, :]
v_samples = full_sample[:, num_repeats:, :]

z_sample_mean = z_samples[:, :, :].mean(axis=1)
z_sample_std = z_samples[:, :, :].std(axis=1)

single_X = cast_array_to_local_type(X_array)
init_mean = f_z_init_mean(single_X)[0]
init_var = f_z_init_var(single_X)[0]

init_vel_var = f_v_init_var(single_X)[0]

print 'Posterior distribution statistics'
print
print 'Initial model: - Mean: ' + str(init_mean)
print '               - Var:  ' + str(init_var)
print
print 'Full HVI model: - Mean: ' + str(z_sample_mean[m.n_hmc_steps])
print '                - Var:  ' + str(z_sample_std[m.n_hmc_steps] ** 2)
print
print 'Velocity model variance: ' + str(init_vel_var)
Posterior distribution statistics

Initial model: - Mean: [ 0.8988  0.5144 -0.0756 -1.8144 -0.6053  1.2234  0.922  -1.0383  0.3988
  1.141   0.0689 -0.2119 -0.9036 -0.6872 -0.0869  0.2268 -1.0339  0.5664
 -1.7155 -1.6149]
               - Var:  [ 0.0266  0.0543  0.1257  0.0449  0.0086  0.0319  0.1534  0.0113  0.193
  0.2356  0.0687  0.0851  0.1472  0.0685  0.0212  0.2111  0.0754  0.1505
  0.0557  0.0521]

Full HVI model: - Mean: [ 1.1074  0.5165 -0.6409 -2.357  -0.6317  1.6201  1.0966 -1.3644  0.4508
  1.2567 -0.2164 -0.1597 -1.121  -1.3046 -0.2349 -0.1019 -1.5526  0.4075
 -1.3688 -2.3848]
                - Var:  [ 0.0519  0.0784  0.1433  0.0518  0.0107  0.067   0.2406  0.0187  0.273
  0.3644  0.0914  0.1393  0.2152  0.1023  0.0435  0.3354  0.1334  0.3178
  0.1508  0.0845]

Velocity model variance: [ 2.5374  0.7248  0.4448  2.4705  2.0529  1.4848  0.3267  1.8692  0.3661
  0.419   0.5219  0.5367  0.4179  0.7956  1.8922  0.3375  0.6573  0.3522
  0.9355  0.964 ]

In [36]:
mean_loss = m.score(repeated_X).mean()
print 'Mean loss for this X: ', mean_loss

single_z_evolution = z_samples[:, 0, :]

R = f_gen_rate(cast_array_to_local_type(single_z_evolution)).astype('float32')

num_steps = m.n_hmc_steps + 1
num_images = num_steps + 1
fig, axs = plt.subplots(1, num_images, figsize=(9*num_images, 9))

img = tile_raster_images(X_array, image_dims, (1, 1), (1, 1))
axs[0].imshow(img, cmap=cm.nipy_spectral)  # cmap=cm.binary

for i in range(num_steps):
    img = tile_raster_images(np.array([R[i]]), image_dims, (1, 1), (1, 1))
    axs[1 + i].imshow(img, cmap=cm.nipy_spectral)
    pot_en_of_image = f_pot_en(single_X, cast_array_to_local_type(np.array([single_z_evolution[i]])))[0]
    axs[1 + i].set_title('Pot. energy: ' + str(pot_en_of_image), fontsize=32)
Mean loss for this X:  117.338691711

Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [39]:
dim1 = 0
dim2 = 1
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [48]:
number_images_per_axis = 10
resolution_per_image = 19  # must be an odd number
x_range = .5
y_range = .5

resolution = number_images_per_axis * resolution_per_image
z_sample_final_mean = z_sample_mean[m.n_hmc_steps]
lower_dim1_limit = z_sample_final_mean[dim1] - 0.5*x_range
upper_dim1_limit = z_sample_final_mean[dim1] + 0.5*x_range
lower_dim2_limit = z_sample_final_mean[dim2] - 0.5*y_range
upper_dim2_limit = z_sample_final_mean[dim2] + 0.5*y_range

latent_array = np.zeros((number_images_per_axis**2, n_latents))

pot_energy_matrix = np.zeros((resolution, resolution), dtype='float32')
x = np.linspace(lower_dim1_limit, upper_dim1_limit, resolution)
y = np.linspace(lower_dim2_limit, upper_dim2_limit, resolution)
for i in range(resolution):
    for j in range(resolution):
        #pos_array = f_z_init_mean(single_X)
        pos_array = np.array([z_sample_final_mean])
        pos_array[0, dim1] = x[i]
        pos_array[0, dim2] = y[j]
        pos_array_of_local_type = cast_array_to_local_type(pos_array)
        pot_energy_matrix[j, i] = f_pot_en(single_X, pos_array_of_local_type)[0]
        if i % resolution_per_image == (resolution_per_image - 1)/2 and j % resolution_per_image == (resolution_per_image - 1)/2:
            latent_array[(i//resolution_per_image) + (number_images_per_axis - 1 - j//resolution_per_image)*number_images_per_axis , :] = pos_array[0, :]

        
print 'Minimum potential energy (at grid points): ' + str(pot_energy_matrix.min())
print 'Maximum potential energy (at grid points): ' + str(pot_energy_matrix.max())

fig, axs = plt.subplots(1, 1, figsize=(15, 15))
CS = axs.contour(x, y, pot_energy_matrix, 40)
plt.clabel(CS, inline=1, fmt='%1.1f', fontsize=10, colors=None)  # colors='black'

F = f_gen_rate(cast_array_to_local_type(latent_array))
img = tile_raster_images(F, image_dims, (number_images_per_axis, number_images_per_axis), (1, 1))

axs.imshow(img, cmap=cm.binary, extent=(x.min(), x.max(), y.min(), y.max()))
axs.set_title('Potential energy surface')
plt.show()
Minimum potential energy (at grid points): 104.446
Maximum potential energy (at grid points): 108.05

In [49]:
resolution = 50
underlying_variance = f_v_init_var(single_X)
velocity_range_for_images = 10.0 * np.sqrt(underlying_variance[0, :])
lower_dim1_limit = np.around(- velocity_range_for_images[dim1])
upper_dim1_limit = np.around(  velocity_range_for_images[dim1])
lower_dim2_limit = np.around(- velocity_range_for_images[dim2])
upper_dim2_limit = np.around(  velocity_range_for_images[dim2])

kin_energy_matrix = np.zeros((resolution, resolution), dtype='float32')
kin_x = np.linspace(lower_dim1_limit, upper_dim1_limit, resolution)
kin_y = np.linspace(lower_dim2_limit, upper_dim2_limit, resolution)
for i in range(resolution):
    for j in range(resolution):
        vel_array = np.zeros((1, m.n_latent)).astype('float32')
        vel_array[0, dim1] = kin_x[i]
        vel_array[0, dim2] = kin_y[j]
        vel_array_of_local_type = cast_array_to_local_type(vel_array)
        kin_energy_matrix[j, i] = f_kin_en(single_X, vel_array_of_local_type)

print 'Minimum kinetic energy (at grid points): ' + str(kin_energy_matrix.min())
print 'Maximum kinetic energy (at grid points): ' + str(kin_energy_matrix.max())

fig, ax = plt.subplots(1, 1, figsize=(9, 9))
CS = ax.contour(kin_x, kin_y, kin_energy_matrix)
#plt.axes().set_aspect('equal', 'datalim')
plt.clabel(CS, inline=1, fmt='%1.1f', fontsize=10)
ax.set_title('Kinetic energy surface')
plt.show()
Minimum kinetic energy (at grid points): 15.9536
Maximum kinetic energy (at grid points): 122.236

Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [50]:
fig, axs = plt.subplots(m.n_hmc_steps + 1, 3, figsize=(18, (m.n_hmc_steps + 1) * 6))
colors = cm.jet(np.linspace(0, 1, 10))

#contour_levels = (198, 200, 202, 204, 206, 208, 210)
#contour_levels = (120, 130, 140, 150, 160, 180, 200, 240, 280)
#contour_levels = (100, 102, 104, 106, 108, 110, 115, 120, 125, 130)
#contour_levels = (400, 402, 404, 406, 408, 410, 412, 416, 420)
#contour_levels = (106, 108, 110, 112, 114, 116, 118, 120, 124, 128)
#contour_levels = (160, 165, 170, 175, 180, 185, 190, 195, 200, 210, 220, 230, 240, 250, 270, 300)
#contour_levels = (174, 175, 176, 177, 178, 180, 182, 184, 186, 190, 200)
#contour_levels = (59, 61, 63, 65, 67, 69, 71, 73, 75, 80, 85, 90)
contour_levels = 40

vel_contour_levels = np.linspace(2.0, 70.0, 18)
#CS0 = axs[0, 0].contourf(x, y, pot_energy_matrix, np.linspace(155, 240, 500))

def colour_for_z_samples(samples):
    mean = samples.mean(axis=0)
    mean1 = mean[dim1]
    mean2 = mean[dim2]
    colour = np.zeros_like(samples[:, 0])
    colour[np.logical_and(samples[:, dim1] < mean1,  samples[:, dim2] < mean2)] = 0
    colour[np.logical_and(samples[:, dim1] < mean1,  samples[:, dim2] >= mean2)] = 2
    colour[np.logical_and(samples[:, dim1] >= mean1, samples[:, dim2] < mean2)] = 4
    colour[np.logical_and(samples[:, dim1] >= mean1, samples[:, dim2] >= mean2)] = 7
    colour[((samples[:, dim1] - mean1) ** 2 + (samples[:, dim2] - mean2) ** 2) < 1e-5] = 9
    return colour.astype('int32')

for i in range(m.n_hmc_steps + 1):
    CS = axs[i, 0].contour(x, y, pot_energy_matrix, contour_levels)
    plt.clabel(CS, inline=1, fmt='%1.0f', fontsize=10)
    axs[i, 0].scatter(z_samples[i,:,dim1], z_samples[i,:,dim2], c=colors[colour_for_z_samples(z_samples[i,:,:])], s=20, alpha=.3, lw=0)
    
    CS_vel = axs[i, 1].contour(kin_x, kin_y, kin_energy_matrix, vel_contour_levels)
    plt.clabel(CS_vel, inline=1, fmt='%1.1f', fontsize=10)
    axs[i, 1].scatter(v_samples[i,:,dim1], v_samples[i,:,dim2], c=colors[colour_for_z_samples(z_samples[i,:,:])], s=20, alpha=.3, lw=0)
    
    pot_energy_distrib = f_pot_en(repeated_X, cast_array_to_local_type(z_samples[i, :, :]))
    if i == 0:
        max_x_value_for_hist = pot_energy_distrib.max() + 5
        min_x_value_for_hist = np.floor(pot_energy_matrix.min()) -5
    pot_energy_distrib_mean = pot_energy_distrib.mean()
    axs[i, 2].hist(pot_energy_distrib, 30, normed=1, range=(min_x_value_for_hist, max_x_value_for_hist))
    axs[i, 2].autoscale(enable=False, axis='both')
    axs[i, 2].axvline(pot_energy_distrib_mean, color='r', linestyle='dashed', linewidth=2)
    axs[i, 2].set_xlim(min_x_value_for_hist, max_x_value_for_hist)
    axs[i, 2].text(pot_energy_distrib_mean + 1.0, 0.8*axs[i, 2].get_ylim()[1], 'Mean: ' + str(pot_energy_distrib_mean))
    axs[i, 1].set_xlim(-velocity_range_for_images[dim1], velocity_range_for_images[dim1])
    axs[i, 1].set_ylim(-velocity_range_for_images[dim2], velocity_range_for_images[dim2])
    axs[i, 1].set_aspect('equal', 'datalim')
    axs[i, 0].set_aspect('equal', 'datalim')

axs[0, 0].scatter(f_z_init_mean(single_X)[0, dim1], f_z_init_mean(single_X)[0, dim2], c='black', s=20)

plt.show()
In [None]:
z_pos = cast_array_to_local_type(np.array([z_sample_mean[m.n_hmc_steps]])) # [z_sample_mean[3]]

second_deriv0 = f_pot_en_2nd_deriv0(single_X, z_pos)
second_deriv1 = f_pot_en_2nd_deriv1(single_X, z_pos)
second_deriv = np.concatenate([second_deriv0, second_deriv1], axis=0)
det = np.linalg.det(second_deriv)
trace = np.trace(second_deriv)

if trace >= 0 and det >= 0:
    eigvalstring = '+ +'
elif det < 0:
    eigvalstring = '+ -'
else:  # so trace < 0 and det >= 0
    eigvalstring = '- -'

print 'Potential energy:'
print f_pot_en(single_X, z_pos)
print
print '1st derivative:'
print f_pot_en_deriv(single_X, z_pos)
print
print '2nd derivative:'
print second_deriv
print
print '2nd deriv. eigenvalue signs: ', eigvalstring
In [None]:
#debugprint(f_z1_by_z0.theano_func.theano_func)
zposss = cast_array_to_local_type(np.array([single_z_evolution[0]]))
print f_z1_by_z0(single_X, zposss, zposss)
print f_zT_by_z0(single_X, zposss, zposss)
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [51]:
np.random.seed(1)

velocity_noise = cast_array_to_local_type(np.random.normal(size=(m.n_hmc_steps, 1, m.n_latent)))
#velocity_noise = np.zeros_like(velocity_noise)
f_z_init_sample
init_pos = f_z_init_sample(single_X) # f_z_init_mean(single_X) # + np.array([0.0, 0.1])
init_vel = f_kin_energy_sample_from_noise(single_X, velocity_noise[0])

num_vels_per_hmc = (m.n_lf_steps + 2)
position_array = np.zeros((m.n_hmc_steps * m.n_lf_steps + 1, m.n_latent))
position_array[0] = init_pos
velocity_array = np.zeros((m.n_hmc_steps * num_vels_per_hmc, m.n_latent))
velocity_array[0] = ma.assert_numpy(init_vel)

for hmc_num in range(m.n_hmc_steps):
    if hmc_num == 0:
        curr_pos = cast_array_to_local_type(init_pos)
        curr_vel = init_vel
    else:
        curr_vel = f_updated_vel_from_noise(single_X, curr_vel, velocity_noise[hmc_num])
        velocity_array[hmc_num * (m.n_lf_steps + 2)] = ma.assert_numpy(curr_vel)
    
    lf_step_results = f_perform_lf_steps(single_X, curr_pos, curr_vel)
    pos_steps = lf_step_results[:m.n_lf_steps]
    vel_half_steps_and_final = lf_step_results[m.n_lf_steps:]
    final_vel = lf_step_results[-1]
    final_pos = pos_steps[-1]
    
    position_array[hmc_num * m.n_lf_steps + 1: (hmc_num + 1)*m.n_lf_steps + 1] = ma.assert_numpy(pos_steps[:, 0, :])
    velocity_array[hmc_num * num_vels_per_hmc + 1: (hmc_num + 1) * num_vels_per_hmc] = ma.assert_numpy(vel_half_steps_and_final[:, 0, :])
    
    curr_pos = final_pos
    curr_vel = final_vel
In [52]:
fig, axs = plt.subplots(1, 2, figsize=(18, 9))
step_color = cm.jet(np.linspace(0, 1, position_array.shape[0]))
CS = axs[0].contour(x, y, pot_energy_matrix, 40)
CS_vel = axs[1].contour(kin_x, kin_y, kin_energy_matrix, vel_contour_levels)
hmc_step_indices = np.arange(0, position_array.shape[0], m.n_lf_steps)
size_array = 40*np.ones((position_array.shape[0],))
size_array[hmc_step_indices] = 100
axs[0].scatter(position_array[:, dim1], position_array[:, dim2], c=step_color, lw=1, s=size_array)
axs[1].set_color_cycle(step_color)

for hmc_num in range(m.n_hmc_steps):
    curr_vel_range = np.arange(num_vels_per_hmc * hmc_num, num_vels_per_hmc * (hmc_num + 1) - 2)
    init_vel_ind = hmc_num * num_vels_per_hmc
    final_vel_ind = (hmc_num + 1) * num_vels_per_hmc - 1
    curr_index = hmc_step_indices[hmc_num]
    next_index = hmc_step_indices[hmc_num + 1]
    for j in curr_vel_range:
        axs[1].plot(velocity_array[j:j+2, dim1], velocity_array[j:j+2, dim2], lw=2)
    axs[1].scatter(velocity_array[init_vel_ind, dim1], velocity_array[init_vel_ind, dim2], c=step_color[curr_index], lw=0, s=100)
    axs[1].scatter(velocity_array[final_vel_ind, dim1], velocity_array[final_vel_ind, dim2], c=step_color[next_index], lw=0, s=100)

for hmc_num in range(m.n_hmc_steps):
    final_vel_ind = (hmc_num + 1) * num_vels_per_hmc - 1
    next_index = hmc_step_indices[hmc_num + 1]
    axs[1].plot(velocity_array[final_vel_ind-1:final_vel_ind+1, dim1], velocity_array[final_vel_ind-1:final_vel_ind+1, dim2], lw=2, c=step_color[next_index])

axs[0].set_aspect('equal', 'datalim')
axs[1].set_aspect('equal', 'datalim')
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [53]:
inv_model_mean_output = np.zeros((m.n_hmc_steps, m.n_latent, num_repeats, m.n_latent))
inv_model_var_output = np.zeros((m.n_hmc_steps, m.n_latent, num_repeats, m.n_latent))

for i in range(m.n_hmc_steps):
    variation_start = z_sample_mean[i] - 2*z_sample_std[i]
    variation_end = z_sample_mean[i] + 2*z_sample_std[i]
    
    time_array = cast_array_to_local_type(np.array([i]).astype('float32'))
    
    for variation_dim in range(m.n_latent):
        z_variation = np.linspace(variation_start[variation_dim], variation_end[variation_dim], num_repeats)
        sample_array = np.tile(z_sample_mean[i], (num_repeats, 1))
        sample_array[:, variation_dim] = z_variation
        local_type_sample_array = cast_array_to_local_type(sample_array)
        if 'velocity' not in m.auxiliary_inv_model_inpt:
            inv_model_mean_output[i, variation_dim] = f_aux_inv_mean(repeated_X, local_type_sample_array, time_array)
            inv_model_var_output[i, variation_dim] = f_aux_inv_var(repeated_X, local_type_sample_array, time_array)
        else:
            velocity_array = 0.0 * local_type_sample_array
            inv_model_mean_output[i, variation_dim] = f_aux_inv_mean(repeated_X, local_type_sample_array, velocity_array, time_array)
            inv_model_var_output[i, variation_dim] = f_aux_inv_var(repeated_X, local_type_sample_array, velocity_array, time_array)
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [54]:
variation_start = z_sample_mean[m.n_hmc_steps] - 2*z_sample_std[m.n_hmc_steps]
variation_end = z_sample_mean[m.n_hmc_steps] + 2*z_sample_std[m.n_hmc_steps]
print variation_start
print variation_end
final_vel_model_mean_output = np.zeros((m.n_latent, num_repeats, m.n_latent))
final_vel_model_var_output = np.zeros((m.n_latent, num_repeats, m.n_latent))

for variation_dim in range(m.n_latent):
    z_variation = np.linspace(variation_start[variation_dim], variation_end[variation_dim], num_repeats)
    sample_array = np.tile(z_sample_final_mean, (num_repeats, 1))
    sample_array[:, variation_dim] = z_variation
    final_vel_model_mean_output[variation_dim] = f_v_final_mean(repeated_X, cast_array_to_local_type(sample_array))
    final_vel_model_var_output[variation_dim] = f_v_final_var(repeated_X, cast_array_to_local_type(sample_array))
[ 0.6519 -0.0435 -1.3981 -2.8122 -0.8383  1.1024  0.1156 -1.6379 -0.5942
  0.0494 -0.8209 -0.9062 -2.0488 -1.9445 -0.652  -1.2602 -2.283  -0.7199
 -2.1454 -2.9661]
[ 1.5629  1.0765  0.1163 -1.9019 -0.4251  2.1377  2.0775 -1.0908  1.4958
  2.4639  0.3881  0.5868 -0.1931 -0.6648  0.1822  1.0563 -0.8223  1.535
 -0.5922 -1.8035]

In [55]:
fig, axs = plt.subplots(m.n_hmc_steps, 2, figsize=(18, 9*m.n_hmc_steps))

for i in range(m.n_hmc_steps - 1):
    axs[i, 0].scatter(inv_model_mean_output[i+1, :, :, dim1], 
           inv_model_mean_output[i+1, :, :, dim2],  
           c=np.transpose(np.tile(np.linspace(0, m.n_latent-1, m.n_latent), (num_repeats, 1))), 
           lw=0, s=5)
    axs[i, 1].scatter(inv_model_var_output[i+1, :, :, dim1], 
           inv_model_var_output[i+1, :, :, dim2],  
           c=np.transpose(np.tile(np.linspace(0, m.n_latent-1, m.n_latent), (num_repeats, 1))), 
           lw=0, s=5)
    axs[i, 0].set_title('Aux. inv. model ' + str(i+1) + ': Mean')
    axs[i, 1].set_title('Final vel. model ' + str(i+1) + ': Variance')
        
if m.n_hmc_steps == 1:
    axs[0].scatter(final_vel_model_mean_output[:, :, dim1], 
               final_vel_model_mean_output[:, :, dim2],  
               c=np.transpose(np.tile(np.linspace(0,m.n_latent-1,m.n_latent), (num_repeats, 1))), 
               lw=0, s=5)
    axs[1].scatter(final_vel_model_var_output[:, :, dim1], 
               final_vel_model_var_output[:, :, dim2],  
               c=np.transpose(np.tile(np.linspace(0,m.n_latent-1,m.n_latent), (num_repeats, 1))), 
               lw=0, s=5)
    axs[0].set_title('Final vel. model: Mean')
    axs[1].set_title('Final vel. model: Variance')
else:
    axs[m.n_hmc_steps-1, 0].scatter(final_vel_model_mean_output[:, :, dim1], 
               final_vel_model_mean_output[:, :, dim2],  
               c=np.transpose(np.tile(np.linspace(0,m.n_latent-1,m.n_latent), (num_repeats, 1))), 
               lw=0, s=5)
    axs[m.n_hmc_steps-1, 1].scatter(final_vel_model_var_output[:, :, dim1], 
               final_vel_model_var_output[:, :, dim2],  
               c=np.transpose(np.tile(np.linspace(0,m.n_latent-1,m.n_latent), (num_repeats, 1))), 
               lw=0, s=5)
    axs[m.n_hmc_steps-1, 0].set_title('Final vel. model: Mean')
    axs[m.n_hmc_steps-1, 1].set_title('Final vel. model: Variance')

plt.show()
In [56]:
print v_samples[0, :, :].mean(axis=0)
print v_samples[1, :, :].mean(axis=0)
print v_samples[2, :, :].mean(axis=0)
print v_samples[3, :, :].mean(axis=0)
print
print v_samples[0, :, :].var(axis=0)
print v_samples[1, :, :].var(axis=0)
print v_samples[2, :, :].var(axis=0)
print v_samples[3, :, :].var(axis=0)
[-0.0539  0.0068 -0.019   0.0661  0.0677  0.0308 -0.0125 -0.1338 -0.0506
 -0.0283  0.0169  0.0002  0.0065 -0.0436 -0.068   0.0012 -0.0338 -0.0032
  0.0205  0.0028]
[ 1.8717 -0.0212 -0.8579 -3.4012 -1.1796  1.6376  0.2313 -1.3531  0.135
  0.2963 -0.2492  0.0159 -0.2276 -1.3966 -0.7814 -0.1425 -0.8697 -0.1587
  0.7156 -2.1059]
[ 0.1567 -0.1525 -0.1131 -1.2434  0.065   0.4919 -0.0324 -0.8482 -0.0952
 -0.088  -0.0084  0.0237 -0.0193 -0.3405 -0.5818 -0.2124 -0.6739  0.0994
  0.3218 -0.8716]
[ 0.2016  0.0198 -0.0798 -1.2251 -0.0706  0.3542  0.0154 -0.4231 -0.0964
  0.0057 -0.1275  0.1449 -0.0653 -0.1648 -0.2539 -0.1611 -0.1546 -0.0976
  0.3707 -0.5758]

[ 2.3796  0.7273  0.4636  2.4914  2.0148  1.4761  0.3233  1.8929  0.3561
  0.4201  0.4935  0.511   0.4069  0.806   1.9455  0.325   0.6227  0.3623
  0.882   0.9238]
[ 2.3386  0.7851  0.5032  3.6804  3.6888  1.5709  0.3254  4.6293  0.3725
  0.3834  0.5878  0.615   0.3528  0.7653  1.6144  0.3023  0.6813  0.2688
  1.0123  0.9103]
[ 2.3836  0.7671  0.4453  2.2305  2.3772  1.3976  0.3394  2.1302  0.3844
  0.3966  0.584   0.603   0.4055  0.801   1.7242  0.2728  0.7011  0.3211
  0.8276  0.9175]
[ 2.5047  0.7706  0.4609  2.6814  2.1707  1.6366  0.3432  1.9547  0.3537
  0.4067  0.5786  0.5281  0.4202  0.761   2.109   0.3128  0.5938  0.3132
  0.931   0.977 ]

In [None]:
final_z_samples = cast_array_to_local_type(z_samples[m.n_hmc_steps, :, :])
final_v_samples = cast_array_to_local_type(v_samples[m.n_hmc_steps, :, :])
final_vel_mean = f_v_final_mean(repeated_X, final_z_samples)
final_vel_var = f_v_final_var(repeated_X, final_z_samples)
final_vel_nll = f_v_final_model_nll(repeated_X, final_z_samples, final_v_samples)
In [None]:
print f_kin_energy_nll(single_X).sum(-1)

print final_vel_nll.mean()
print final_vel_nll.min()
print final_vel_nll.max()
In [None]:
#m.parameters.view(m.auxiliary_inv_model.mlp.layers[2].weights)
print m.parameters._var_to_slice[m.auxiliary_inv_model.mlp.layers[2].weights]
#grad[594190:594990]
#old_weights_array = m.parameters.view(m.auxiliary_inv_model.mlp.layers[0].weights)
#new_array = old_weights_array.copy()
#new_array[-1, :] = 1.0 + old_weights_array[-1, :]
#m.parameters.__setitem__(m.auxiliary_inv_model.mlp.layers[0].weights, new_array)
print m.parameters.view(m.auxiliary_inv_model.mlp.layers[2].weights).shape
In [None]:
np.set_printoptions(suppress=True)
print np.around(np.reshape(mean_grad[396390:553790], (787, 200)), 5)[:198]  # 396390:553790
np.reshape(mean_grad[594190:594990], (200, 4))  # 396390:553790
#print np.around(np.reshape(m.parameters.data[396390:553790].astype('float32'), (787, 200)), 5)[:50]
Node Commands Syntax: node {operator} [options] [arguments] Parameters: /? or /help - Display this help message. list - List nodes or node history or the cluster listcores - List cores on the cluster view - View properties of a node online - Set nodes or node to online state offline - Set one or more nodes to the offline state For more information about HPC command-line tools, see http://go.microsoft.com/fwlink/?LinkId=120724.
In [None]:
fig, axs = plt.subplots(4, 2, figsize=(18, 36))
# TODO: Analysis of how final_vel_mean and final_vel_var depend on z (since they all share the same x)

print z_samples[3, :, :].mean(axis=0)
print z_samples[3, :, :].var(axis=0)
print v_samples[3, :, :].mean(axis=0)
print v_samples[3, :, :].var(axis=0)
print f_v_init_var(np.array([X[X_index, :]]))

print final_vel_nll.mean()
plt.boxplot(final_vel_nll, whis=1)
plt.show()
In [None]:
centers = np.zeros((10,n_latents))
stddevs = np.zeros((10,n_latents))
centers_init = np.zeros((10,n_latents))
stddevs_init = np.zeros((10,n_latents))
for i in range(10):
    Li = f_z_sample(X[Z.argmax(1) == i])
    centers[i] = Li.mean(axis=0)
    stddevs[i] = np.std(Li, axis=0)
    
    Li_init = f_z_init_sample(X[Z.argmax(1) == i])
    centers_init[i] = Li_init.mean(axis=0)
    stddevs_init[i] = np.std(Li_init, axis=0)
In [None]:
fig, axs = plt.subplots(1, 2, figsize=(18, 9))
axs[0].scatter(centers[:, dim1], centers[:, dim2], c=range(10), s=50)
axs[0].scatter(centers_init[:, dim1], centers_init[:, dim2], c=range(10), s=50, marker=u's')

axs[1].scatter(centers[:, dim1], centers[:, dim2], c=range(10), s=50)
axs[1].scatter(centers[:, dim1] + stddevs[:, dim1], centers[:, dim2], c=range(10), s=50, marker=u'>')
axs[1].scatter(centers[:, dim1] - stddevs[:, dim1], centers[:, dim2], c=range(10), s=50, marker=u'<')
axs[1].scatter(centers[:, dim1], centers[:, dim2] + stddevs[:, dim2], c=range(10), s=50, marker=u'^')
axs[1].scatter(centers[:, dim1], centers[:, dim2] - stddevs[:, dim2], c=range(10), s=50, marker=u'v')

#axs[0].set_xlim(-1.2, 1.2)
#axs[0].set_ylim(-1.2, 1.2)
#axs[1].set_xlim(-1.2, 1.2)
#axs[1].set_ylim(-1.2, 1.2)

print (centers[:, dim1] - centers_init[:, dim1])
print (centers[:, dim2] - centers_init[:, dim2])
print (stddevs[:, dim1] - stddevs_init[:, dim1])
print (stddevs[:, dim2] - stddevs_init[:, dim2])
In [None]: